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CHAPTER  I 
INTRODUCTION 

When  a  practical  problem  in  science  or  technology  permits 
mathematical  formulation,  the  chances  are  rather  good  that  it 
leads  to  one  or  more  differential  equations.   This  is  certainly 
true  of  the  vast  catagory  of  problems  associated  with  beams  and 
columns,  heat  transfer  and  fluid  flow,  elasticity  and  electricity, 
etc. 

Many  of  these  problems  which  can  be  neatly  formulated  as 
differential  equations  can  go  no  further  for  lack  of  solutions. 
At  this  point  one  discovers  how  few,  relatively  speaking,  are  the 
equations  that  have  solutions  in  closed  form. 

Faced  with  this  situation,  various  numerical  methods  have 
been  devised  that  squeeze  the  desired  information  out  of  the 
differential  equation  directly.   However,  most  of  these  methods 
yield  a  solution  over  a  limited  range  of  the  domain,  and  this  is 
in  the  form  of  a  table,  giving  values  of  the  dependent  variable 
and  sometimes  its  few  derivatives  for  specific  values  of  the  in- 
dependent variable. 

Many  times  this  defines  the  solution  well  enough,  particu- 
larly if  the  solution  has  been  tabulated  with  a  sufficiently 
small  increment  of  the  independent  variable.   But  however  fine 
the  increment  be,  the  solution  does  not  compare  with  an  analytical 


one,  in  that  it  cannot  be  handled  analytically  -  put  in  equations, 
differentiated  or  integrated  as  desired,  and  so  on. 

Iterative  methods,  Improving  an  assumed  analytical  solution 
with  each  step,  have  been  suggested,  and  successfully  applied  to 
problems  of  initial-value  type  and  also  to  those  boundary  value 
problems  for  which  Green's  function  is  known.   Iteration  was  first 
applied  to  a  technical  eigenvalue  problem  in  I898  by  L.  Vianello 
[21]*  in  a  study  of  buckling  problems.   The  process  was  applied 
by  A.  Stodola  [18]  in  190^  to  the  problem  of  critical  speeds  of 
rotating  shafts.   The  theory  of  the  method  had  already  been  pre- 
sented by  H.  A.  Schwarz  [16]  in  1885  and  been  developed  by  E. 
Picard  [13]  .   It  is  interesting  to  note  that  all  these  early  works 
dealt  with  a  particular  class  of  eigenvalue  problems.   No  effort 
was  ever  made  to  apply  the  technique  to  equilibrium  problems  and 
not  many  instances  could  be  cited  when  an  attempt  was  made  to 
develop  this  technique  further  to  apply  to  a  wider  class  of  pro- 
blems. 

It  has  been  attempted  in  the  present  work  to  develop  a 
technique  for  obtaining  an  analytical  solution  to  boundary  value 
problems  governed  by  quasi-linear  differential  equations.   The 
solution  was  intended  to  be  a  power  series  in  the  independent 
variable.   The  treatment  here  was  limited  to  problems  in  one 
dimension,  which  include  the  initial  value  problems  as  a  sub- 
class.  However,  the  author  has  a  feeling  that  with  enough  modi- 
fications, the  technique  could  be  extended  to  problems  in  more 
than  one  dimension. 


*  Numbers  in  brackets  refer  to  references  in  Bibliography 


CHAPTER  II 
THE  PROBLEM 
Nomenclature 

English  alphabets: 

A  area  of  cross-section 

b  breadth  of  section 

B  '  coefficient  matrix  in  the  matrix  equation  representing  linear 

boundary  conditions 
Bi  ordinary  differential  operator 

c,  coefficients  in  the  series  expansion  of  the  solution  function 
C  column  vector  of  c^>  i  =  1»  2,...,  m. 

d  depth  of  section,  dr.  <*2  ""  end  depths 

D  domain  of  definition;  d/dx;  coefficient  matrix  in  DC  =  F 

(refer  p.  13) 
E  modulus  of  elasticity 

F  constant  column  vector  in  DC  =  F  (refer  p.  13) 
G  ordinary  differential  operator 
I  moment  of  area  of  cross-section 
K  constant  column  vector  in  the  matrix  representation  of  linear 

boundary  conditions 
L  length 

m  order  of  the  governing  differential  equation 
M  quasi-linear  ordinary  differential  operator,  containing  the 

highest  order  derivative  term 
n  number  of  boundary  conditions  at  the  left  end 


N  ordinary  differential  operator 
x  the  independent  variable 
y  the  dependent  variable 

Greek  alphabets: 

1/4 

a  ratio  of  end  depths  =  d2/d1;  X 

3  X^VL 

Y  (cos  a  -  cosha)/(  sinha  -  sin  a) 

X   eigenvalue,    in  the  non-dimensional   formulation 

i 

p   density 

I   summation  sign 

id  angular  frequency 

Subscripts  and  superscripts: 

Number! cal  subscript  to  y  or  x  represent  that  particular  mode. 
Roman  superscript  to  y  represents  differentiation. 
Numberical  or  alphabetic  superscript  in  parentheses  represent  a 
particular  iteration. 


Statement  of  the  Problem 

A  boundary  value  problem,  in  general,  may  be  stated  as 

M[y]  =XN[y]         in  D 

(II-D 
B*[y]    =   0  i  =  1,  2,.  ..  ,m 

on  the  boundary  of  D. 

Here  D  is  the  domain  of  definition  of  the  problem  and,  in 
case  of  two-point  boundary  value  problems,  is  a  one-dimensional 
continum.   M  and  N  are  ordinary  differential  operators,  M  being 
quasi-linear  in  nature.   The  order  of  the  differential  equation 
m  is  the  same  as  that  of  the  operator  M  and  is  larger  than  the 
order  of  the  operator  N. 

The  presence  or  absence  of  X  ,  an  undetermined  parameter, 
determines  whether  the  problem  belongs  to  the  eigenvalue  class  or 
equilibrium  class.   In  the  case  of  eigenvalue  problems  the  boun- 
dary conditions  are  essentially  homogeneous.   The  formulation  is 
satisfied  by  an  infinity  of  values  of  X  and  the  corresponding 
eigenf unctions  y. 

In  the  present  work  a  solution  of  the  formulation  was  sought 
in  the  shape  of  y  expressed  as  power  series  in  the  independent 
variable  x.   In  case  of  eigenvalue  problems,  the  first  few  eigen- 
values and  the  corresponding  mode  shapes  constituted  what  was 
expected  as  solution. 


The  Technique 


The  governing  differential  equation  of  the  formulation, 
being  quasi-linear  in  nature,  was  very  easily  solved  for  the 
highest  derivative  of  y. 


Dmy  =  Gly] 

or  y  =  D-m(G[y])  ♦  Jo.  x1'1  (H-2) 

i=l 

where  D  =  d/dx  and  the  negative  powers  indicate  integration.  The 
c.,  the  constants  of  integration,  were  determined  by  the  use  of 
the  m  boundary  conditions  - 

B.  [y]  =  0  i  =  1,  2 m. 

An  iterative  process  was  then  set  up 


m 

I 
i=l 

(r+1),  _n  (II"3) 


B1  [yVJ-  iy]  =  0 

I*   =   Xf    b|  *  M  1 

To  start  with,  a  polynomial  was  selected  as  an  initial  guess 
for  y  i.e.  y^  .   Any  polynomial  would  serve  the  purpose  but  it 
is  advisable  to  select  one  that  satisfies  all  the  boundary  con- 
ditions.  The  proper  choice  of  y^  '  would  certainly  accelerate 
the  convergence  of  the  process. 

This  was  enough  to  proceed,  if  the  parameter  X   did  not  enter 
the  formulation.   Otherwise,  it  sometimes  necessiated  a  guess  for 


X  (in  addition  to  that  of  solution  function)  at  every  step  of 
iteration.   However,  this  was  easily  furnished  by  Raleigh's  quo- 
tient of  the  trial  function. 

It  was  observed  that  the  iteration  process  here  built  up  a 
power  series  approxitation  to  the  true  solution,  each  iteration 
adding  substantially  a  few  terms  to  the  expansion  (usually  not 
less  than  the  order  of  the  differential  equation,  m.).   If  the 
process  was  to  converge  (as  it  did  and  could  be  expected  to  in 
most  of  the  well  behaved  cases)  after  enough  iterations  the  first 
few  terms  would  quit  changing.   The  convergence  could  be  observed 
by  comparing  the  successive  iterates. 

In  eigenvalue  problems  the  process  converged  to  the  mode 
corresponding  to  the  smallest  eigenvalue.   The  orthogonality  con- 
dition was  used  to  extract  the  higher  modes. 
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Convergence  of  the  Iteration 

For  initial-value  problems,  the  scheme  is  the  same  as  the 
extension  of  Picard's  method  for  higher  order  equations.   The 
proof  has  been  exhibited  in  the  references  [15]  »  [12] .   Levy  and 
Baggott  [8]  give  a  proof  of  convergence  for  systems  with  second 
order  differential  equations  of  both  initial-value  and  equilibrium 
type.   A  general  proof  to  cover  eigenvalue  problems  cannot  be 
given.   However,  a  proof  for  convergence  of  an  iteration  proce- 
dure in  general  is  given  by  Collatz  (  [2],  pp.  36-^-8).   Also,  the 
convergence  of  an  iteration  procedure,  employing  the  inverse  of 
the  operator  M,  is  suggested  in  reference  (  [3],  p.  30  3)-   Vy'hen  M 
is  simply  Dm»  the  method  presented  in  the  thesis  coincides  with 
the  method  mentioned  above.   The  proof  is  given  below: 

The  eigenvalue  problem 


M[y]  =  XN[y] 


(II-*) 


B1[y]  =0,       1=1,2 m 


can  be  written  as 


y  =*G[y] 

where  G  =  M  N  . 

Assume  that  the  expansion  theorem  holds,  and  also  that 
A,  <  \     <_   \„  <_   ju  <,  •••   Any  admissible  function  y*    can  then 

be  expressed  as  a  linear  combination  of  the  individual  modes  y^. 

(0)  r 

y    =  clyl  +  I    ciyi 
i=2 


Starting  from  an  arbitrary  admissible  function  y^   ,  an 
iteration  process  is  set  up  according  to  the  recurrence  relation 

y(n+l)  =G[y(n)] 

n  =  1,  2,  3, - 

For  an  eigenfunction  y, 

y(1)    -W0)l 


=  G[ciyi  +  Ivi 


i=2 

c, 


=  x     yi  +   I    l1  yi 
Al     -1        i=2Ai     x 

oo         X 


y<2>    ^Gly'1') 


so       X 

=  G[f(cy     +    I    ^.y.M 
Al     -1  -1        i=2Ai  1   1 

=  t  (-r  yi  +  l  -r  -r  yi> 

Al      Al     x        i=2Ai   Ai      x 

-t4<Vi*lfa<^'i) 
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Continuing  in  the  same  manner, 

y<n)  •■h^m*  V&\'t>- 


i=2  1 

**  n 
As  n  gets  larger,  the  quotient  (—)   tends  to  zero.   So  that, 

\ 


(n)  _  .1 

(II-5) 


y    =7-Hclyl 
Al 


and  Xt  ■  - 


(n) 


1  "  y(n+l)  • 

Thus  the  iteration  converges  to  the  eigenfunction  corres- 
ponding to  the  smallest  eigenvalue.   In  general,  the  process  will 
converge  to  the  eigenfunction  corresponding  to  the  first  non-zero 
ci,  if  carried  out  exactly.   While  working  approximately,  however, 
small  components  of  y-  will  be  inevitably  introduced,  and  the 
process  will  finally  converge  to  the  eigenfunction  corresponding 
to  the  smallest  eigenvalue  independent  of  the  choice  of  y^  ' . 


CHAPTER  III 
COMPUTERIZATION 

It  is  difficult  to  estimate  in  advance  how  much  computation 
will  be  required  to  obtain  a  solution  by  iteration.   The  amount 
of  computation  per  cycle  increases  with  the  number  of  cycles. 
The  number  of  cycles  required  depends  on  the  accuracy  desired  and 
on  the  particular  system  to  be  solved.   So  that,  at  some  stage, 
hand  computation  goes  out  of  the  question  and  recourse  has  to  be 
taken  to  some  high-speed  computing  machine.   The  problems  illus- 
trated in  the  present  work  were  programmed  for  IBM  1 410/7010 
Operating  System  (1410-PR-155)  F0RTRAN-1410-F0-970. 

Standardization 

For  adapting  the  technique  to  digital  computers,  standardi- 
zation of  the  formulation  was  done.   First  the  formulation  was 
reduced  to  non-dimensional  form.   One  of  the  two  boundary  points 
was  made  to  coincide  with  the  zero  of  the  independent  variable 
axis  ( x-axis) ,  and  the  scale  of  the  variable  x  was  adjusted  so 
that  the  separation  between  the  boundaries  was  unity,  unless 
otherwise  necessary. 
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Boundary  Conditions 

Each  iteration  involved  a  m-fold  integration  and  hence  gave 
rise  to  m  constants  of  integration.   These  constants  were  to  be 
determined  so  as  to  make  the  iterate  satisfy  the  boundary  condi- 
tions.  Application  of  each  condition  resulted  in  an  algebraic 
relation  in  the  m  constants  of  integration  to  be  determined. 
Thus  a  set  of  m  algebraic  equations  was  obtained.   The  complexity 
of  these  equations  naturally  depended  upon  that  of  the  particular 
boundary  conditions.   Fortunately  for  the  computer,  all  the  pro- 
blems used  for  illustration  happened  to  possess  linear  boundary 
conditions.   Nevertheless,  this  was  not  a  limitation  to  the  use- 
fullness  of  the  scheme,  but  it  did  reduce  the  amount  of  computa- 
tion required  for  fixing  the  m  constants.   Also,  most  problems 
of  practical  importance  have  linear  boundary  conditions. 

The  most  general  linear  boundary  conditions  could  be  repre- 
sented by  the  matrix  equation 

B  I  =  K 

where  B  is  the  coefficient  matrix,  Y  is  the  column  matrix  with  y, 
y  »  y   »  y   t .  y      for  its  elements,  and  K  is  the  constant 


column. 


Nc 

y  =  I  c.x1-1 

i=l  1 


i=p+l  Cl-P-1)!   i 
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By  this  time  all  but  the  first  m  coefficients  c±   were  known. 
Hence,  at  boundary  points,  where  x  is  specified,  y  and  any  of  its 

derivatives  were  just  some  linear  combination  of  the  ci (1=1,2, , 

m) .   Hence,  the  unknown  column  Y  could  be  written  as 

Y  =  S  C  +  T 

where  C  is  the  column  of  undetermined  coefficients  cit  and  S  and 
T  were  respectively  the  square  and  column  matrices  determined  as: 


s« 


0  ,  J  <  i 


t  =  r   (,1-1);  x>i-l 


After  substitution 

B  (SC  +  T)  =  K 

BSC  =  K  -  BT 

and  C  =  (BS)"1  (K  -  B  T) 
=  DF 

where  D  =  (BS)"1  (III-D 

P  =  K  -  BT. 

The  subroutine  DMTRX  was  programmed  to  find  the  D  matrix 
while  subroutine  BOUND  constructed  the  column  vector  F.   Finally, 
the  constants  c,  were  obtained  by  matrix  multiplication. 


Ik 


The  problem  of  obtaining  corresponding  matrices  for  obtaining 
modes  corresponding  to  higher  eigenvalues  is  discussed  later 
(Intermediate  Eigenvalues  p.  34) 
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Subroutines 

A  variety  of  operations  were  expected  to  be  performed  re- 
peatedly.  Subroutines  were  programmed  so  that  they  could  be  in- 
serted in  any  other  programs  and  called  for  whenever  desired. 
These  subroutines  along  with  their  limitations  are  briefly  dis- 
cussed below. 

i)  Subroutine  DIFFER  (C,NC,D,ND,N) :  This  subroutine  was  used  to 
differentiate  a  given  polynomial  C  having  No  terms,  N  times  to 
obtain  a  polynomial  D.  Nd,  the  number  of  terms  in  D  equaled  Nc 
minus  N.   The  only  restriction  on  N  was  that  it  should  be  non- 
zero positive  integer. 

ii)  Subroutine  INTGRA  (C,NC,D,ND,N) :  This  subroutine  gave  a  poly- 
nomial D  as  the  Nth  integration  of  the  polynomial  C.   Nd  equalled 
Nc  plus  N.   N  should  be  a  non-zero  positive  integer.   The  con- 
stants of  integrations  would  be  left  undetermined. 
iii)  Subroutine  SSMUL  (A,NA,B,NB,C,NC) :  The  polynomials  A  and  3 
were  multiplied  to  get  polynomial  C  by  this  subroutine.   Na  and 
Nb  both  should  be  non-zero  positive  integers. 

iv)  Subroutine  SEDIV  (A,NA,B,NB,C,NC) :  This  subroutine  was  meant 
to  obtain  the  polynomial  C  as  the  division  of  polynomial  A  by 
polynomial  B.   Nc,  the  number  of  terms  in  C,  should  be  specified 
and  the  process  is  terminated  after  calculating  that  many  terms. 
v)  Subroutine  INVRS  (A,N):  This  subroutine  was  designed  to  find 
the  inverse  of  a  square  matrix  A  of  order  N  by  Gauss-Sidel  re- 
duction process.   In  case  the  inverse  did  not  exist  the  program 
had  instruction  to  write  out  'NO  INVERSE1. 
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vi)  Subroutine  MATMUL  (A.N.3) :  This  subroutine  was  specifically 
made  to  multiply  a  square  matrix  A,  of  order  N,  by  a  column  matrix 

B. 

vii)  Subroutine  INTGRD:  The  aim  of  the  subroutine  INTGRD  was  to 
construct  G[y^r'l  from  the  rth  iterate.   Tne  sequence  of  in- 
structions solely  depended  upon  the  operator  M  (and  N,  if  pre- 
sent) of  the  specific  problem.   This  subroutine,  in  case  of 
eigenvalue  problems,  normalized  the  iterate  every  time  and  also 
computed  the  Raleigh's  quotient,  if  necessary. 

viii)  Subroutine  INITL:  This  subroutine  inilized  the  particular 
problem. 
ix)  The  main  program:  This  program  controlled  the  order  in  which 

all  other  subroutines  were  called. 

x)  Subroutine  RESULT:  This  subroutine  writes  out  what  is  supposed 

to  be  the  answer. 


CHAPTER  IV 
ILLUSTRATIVE  PROBLEMS 

1)  Equilibrium  Problems: 

i)  A  Beam  on  Elastic  Foundation:-  The  equilibrium  of  a  flexible 
beam  subjected  to  a  uniformly  distributed  load,  while  resting  on 
a  continuous  elastic  foundation  was  considered.   The  non- 
dimensional  formulation  as  taken  from  'Engineering  Analysis'  by 
Crandall  (  [3],  p.  195)  was 

yiv  +  y  =  1  0  <  z^l 


and  y(Q)  =  y  (Q) 


=  0  (IV-1) 


y(l)  =y  (1)  =  °   * 

In  problems  of  this  type,  which  possess  unique  non-zero 
solutions,  the  iteration  could  be  started  from  an  initial  guess 

,<9>  .  6. 

In  the  reference  the  problem  was  solved  in  several  ways. 
The  solution  by  Ritz's  stationary  functional  method  with  trial 
family  (  [3],  p.  235) 

y  =  c1  x  (1-x)  +  c2  x2  (1-x2) 
was  taken  for  comparison.   The  solution  contained  only  five  terms, 
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did  not  satisfy  all  the  boundary  conditions,  and  therefore  is  not 
exact.   It  is  compared  with  the  first,  third  and  sixth  iterates, 
in  Table  I. 


Table   I 


Results:    A  Beam  on  Elastic   Foundation 
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Coeff. 

Number 

Comparison 
Solution 

First 

Iterations 

Third 

Sixth 

1 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

2 

0.041249 

0.041667 

0.041249 

0.041249 

3 

0.000032 

0.000000 

0.000000 

0.000000 

4 

-0.082434 

-0.083333 

-0.082646 

-0.082646 

5 

0.041217 

0.041667 

0.041667 

0.041667 

6 

-0.000344 

-0.000344 

7 

0.000000 

0.000000 

8 

0.000098 

0.000098 

9 

-0.000025 

-0.000025 

10 

0.0000001 

0.0000001 
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ii)  Bending  of  a  Strut:-  This  problem  with  its  non-dimensional 
formulation  was  taken  from  "The  Numerical  Treatment  of  Differen- 
tial Equations"  by  Collatz  (p.  143) .   It  considered  the  bending 
of  a  strut  with  varying  flexural  rigidity,  and  axial  compressive 
load,  by  a  distributed  transverse  load.   The  equations  of  equili- 
brium were 

y11  +  (1  +  x2)  y  +  1  =  0    -I  <  i  <  1 

(IV-2) 

y<-u  =  y(+D  =  °- 

To  reduce  this  formulation  to  the  desired  form,  a  change  of 
variable  was  made. 

2  x  =  x  +  1.  (IV-3) 

With  this 

y11  +  4  (1  +  (2  -  4  x2)  y)  =  0    0  <  x  <_  1 

(IV-4) 

*<0)  =  y(D  =  °' 

The  first  seven  terms  of  the  power  series  solution  of  the 
formulation  (IV-2)  were  available  in  the  reference  (p.  225) .   For 
comparison,  the  change  of  variable  (IV-3)  was  made.   The  solution 
appears  in  the  second  column  of  Table  II. 

It  may  be  noted  that  the  initial  guess 


(0)  m 
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did  not  satisfy  any  of  the  boundary  conditions.   This  did  not  in 
anyway  prevent  convergence,  since  at  the  very  next  opportunity 
the  boundary  conditions  were  forced  on  the  iterate.   The  results 
are  given  in  Table  II. 
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Table  II 


Results:  Bending  of  a  Strut 


Comparison 

Coef f . 

Reference 

Iterations 

Number 

Solution 

Second 

Seventh 

Tv/elfth 

1 

0.775  x  10~8 

0.00000 

0.00000 

0.00000 

2 

3.47293 

4.66667 

3.48506 

3.47318 

3 

-1.99998 

-6.00000 

-2.00000 

-2.00000 

4 

-4.63078 

2.66667 

-4.66578 

-4.63129 

5 

+5.96531 

-1.33333 

5.99912 

5.96463 

6 

-2.53309 

-2.51655 

-2.52592 

7 

-2.96677 

-3.05473 

-2.99460 

8 

4.43757 

4.57156 

4.51827 

9 

-1.80832 

-1.99472 

-1.99803 

10 

-0.96844 

-0.65345 

-0.60670 

11 

2.068897 

1.55658 

1.51371 

12 

-1.58388 

-0.91128 

-0.90380 

13 

0.731512 

0.05860 

0.07669 

14 

-0.21578 

0.27760 

0.26406 

15 

0.03083 

-0.20495 

-0.21599 
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2) 

Eigenvalue  Problems: 

i) 

Vibrations  of  a  Beam 

Fixed 

at  One  End  and  Hinge 

d  at 

Another: - 

This  problem,  being  ver;y 

•  common,  could  be  found  in 

almc 

st  every 

text  on  vibration  of  elastic 

bodies.   The  present 

formulation 

appeared  on  page  255  of 

the  reference  [  20]  . 

Fiv  =  3^ 

0  1  x  1  L 

y(o)  =  yl(o)  = 

=  0 

(iv-5) 

y(D  =  y(L)  = 

0 

To  make  it  dimensionless 

,  new  variables  were 

introduced- 

y  =  y/L 

(IV-6) 

x  =  x/L. 

The  new  formulation  was 

yiv  =  *  y 

0  <  x  <  1 

y(o)  =  y(o)  = 

0 

(IV-7) 

y(D  =  y(i)  = 

0. 

and  X  =  3T, 

The  exact  solution 

to  this  formulation  is 

=  (cos  ax  - 

cosh 

ax)  +  y(sln  <*x  -  sinh  ax) 
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i 

4  _ 


where  a  =  A 

cos  a  -  cosh  a 


Y> 


sinh  a  -  sin  a 


The  a i s  are  given  by  the  equation 
Tan  a  =  Tanh  a  , 


The  exact  solution  corresponding  to  the  smallest  eigenvalue 
was  expanded  in  a  power  series  and  the  various  iterates  compared 
with  it.  (Table  III) . 
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Table  III 

Results : 

Beam  Vibrations 

-  First  Mode 

Coeff. 
Number 

Comparison 
Solution 

Iterations 
Tenth         Fifth 

Second 

3 

1.0000000 

1.0000000 

1.0000000 

■ 
1.0000000 

4 

-1.3098845 

-1.3098847 

-1.3098874 

-1.3157893 

7 

0.6603361 

0.6603360 

0.6603518 

0.7368418 

8 

-0.3706989 

-0.3706989 

-0.3707157 

-0.5263155 

9 

0.0000000 

0.0000000 

0.1052631 

11 

0.0311460 

0.0311460 

0.0311551 

12 

-0.0111266 

-0.0111266 

-0.0111328 

15 

+0.0003082 

0.0003082 

0.0003095 

16 

-0.0000807 

-0.0000807 

-0.0000814 

19 

0.0000010 

0.0000010 

0.0000011 

20 

-0.0000002 

-0.0000002 

-0.0000003 

h 

0.00420661 

0.00420651 

0.00376984 

All 

(4n+l)st  and  (4n+2)nd  coeff. 

are  zeros  except  if  last 

term  of  iteration. 

Comp 

•arison  solution 

is  series  expansion  of 

cosS, -coshs. 

(cosSTX-coshSnX)  +  — r-r-r .     (sinB,x-s 

1       1     sinnX, -sinX,       1 

linns,  x) 
1 

with 

i  61  =  3.9266023 

■  X.  ,  the  ei 

.genvalue  obtained  from  the 

tenth  iteration. 
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ii)  Transverse  Vibrations  of  a  Tapered  Beam  of  Rectangular  Cross- 
section,  Simply  Supported:-  This  problem  was  studied  by  the 
author  as  an  assignment  for  the  course  'Machine  Vibrations  II*. 
The  derivation  of  the  formulation  follows: 

In  addition  to  the  usual  assumption  of  small  deflections, 
it  was  assumed  that  the  taper  is  very  small:  so  that  the  differ- 
ential equation  of  motion  can  be  approximated  by  that  for  a  uni- 
form beam  (see  formulation  (IV-5))« 

yiY  =  By         0  <_  x  <_  L 

with  end  conditions 

y<c»  ="('o)  -° 

Onlytthe  term  3  =  &sw   ([20],  p.  255)  is  no  longer  constant 

since  A  and  I  are  variable  for  the  problem  considered  here,  b, 
width  of  the  section  was  assumed  to  be  uniform  and  the  depth  was 
assumed  to  vary  linearly  with  x. 

d(x)  =  dl  +  (d2-di}i      • 

=  d1[l  +  (a-l)-Xr] 


where  d,  and  dp  are  depths  at  left  and  right  ends  and  a 


d2/dl« 
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A(x) 

_  b  d(2) 

b  d?  v/12 

12 

d12[l+(a-l)|]2 

34  = 

2 
12pu       1 

Edx2   l+(a-l)£2 

With  this 

the  differential 

equation  becomes 

l+(a- 

2 

- . x  2  — iv    12 pu   — 

■1}L    y    "     2   Y- 

u                    Ed^ 

With  the  change  of  variable  (IV-6) 

[l+(a-l)x]ylv  =  Xy 

0  <  x  <_  1 

y(0) 

"  *<o)  "  ° 

(IV-8) 

y(D 

=  y(D  "  ° • 

- 

where  X  -  12q^ 
Ed,1* 

The  finite  difference  method  with  ten  intervals  was 

used  to 

solve  the  problem.   For  various 

values  of  a ,  the 

modes  correspond- 

ing  to  the  smallest  eigenvalues 

were  determined. 

These  solutions 

formed  the  basis  for  comparison 

for  the  answers 

given  by 

the  pre- 

sent  technique. 

• 
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Good  agreement  was  found  for  solutions  with  values  of  a  = 

1,  1.1, ,1.9.   But  for  a   =  2  the  process  apparently  converged 

to  a  solution  which  did  not  agree  with  that  obtained  by  the 
finite  difference  method.   An  investigation  into  the  matter  re- 
vealed that  the  reason  was  in  the  divergence  of  the  power  series 
expansion  of  l/(l+x)  .   Forming  this  series  was  an  important  step 
in  constructing  the  integrand  at  each  iteration.   As  might  be  ex- 
pected the  error  was  corrected  by  turning  the  beam  about  its 
centre.   The  value  of  a   was  now  1/2  and  the  expansion  of  1/(1-0. 5x) 
converges  rapidly  enough.   The  various  solutions  are  listed  in 
Table  IV,  along  with  the  solutions  by  the  finite  difference  method. 

It  was  expected  that  the  analytical  solution  would  be  more 
accurate,  since  the  difference  method  used  only  ten  intervals,  and 
inversion  of  a  matrix  of  order  nine  was  involved.   This  was  easily 
verified  for  a  =  1,  (see  Table  V)  that  is  a  uniform  beam.   For 
this  case  the  mathematically  exact  solution  was  known.   It  was 
found  that  the  answers  from  the  difference  method  were  a  bit  too 
low.   There  was  no  way  to  compare  the  accuracies  of  the  two 
methods  in  case  of  other  values  of  a.   But  one  thing  was  certain, 
that  the  same  trend  (the  solutions  by  the  present  technique  being 
slightly  larger  than  those  by  the  difference  method)  was  followed 
throughout. 


Table  IV 


Results:  Vibrations  of  Tapered  Beam 
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Dia 


If 


Ratio  Difference  method  lp 

a      Comparison  Solution   Polynomial  exp, 


(XlP"Xlf)lQQ 
A4p 


1.0 

3.1280996 

3.1415925 

1.1 

3.2051857 

3.2183947 

1.2 

3.2785025 

3.2920185 

1.3 

3.3490291 

3.3628388 

1.4 

3.4170429 

3.4311626 

1.5 

3.4828436 

3.4972445 

1.6 

3.5466312 

3.5612992 

1.7 

3.6085762 

3.6235085 

1.8 

3.6688213 

3.6835085 

1.9 

3.7275184 

3.7123093 

2.0 

3.7847778 

*3. 8005053 

0.41039377 
0.41042200 
0.41056877 
0.41065602 
0.41151387 
0.41177847 
0.41187216 
0.41188344 
0.39872854 
-0.40969377 
0.41382655 


*This  value  was  obtained  by  turning  the  beam  about  its  center 
and  working  the  problem  with  a  =   0.5. 


Table  V 

Results:    Tapered  Beam  Problem  with 
a  =  1      i.e.    a  uniform  beam. 
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Coeff . 

Number 


Exact  Solution 
Comparison  Solution 


Present  Solution 
Tenth  Iteration 


Al 
2 

4 

6 

8 

10 

12 

14 

16 

18 

20 


1.0000000 

-1.6494792 

0.8162345 

-0.19233?^ 
0.0264380 

-0.0023787 
+0.0001509 
-0.0000071 
0.00000025 
Less   than  10 


-8 


3.1415925 

1.0000000 

-1.6449338 

+0.8117423 

-0.1907517 

+0.0261478 

-0.0023460 

0.0001484 

0.0000070 

0.00000025 

Less  than  10 


-8 


Note:-  All  odd  coefficients  in  both  solutions  are   zeros. 
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ill)  Longitudinal  Vibrations  of  a  Cantilever  of  Varying  Cross- 
Section:-  This  eigenvalue  problem  was  formulated  and  solved  in 
reference  (  [2  ] ,  p.  1^7).   The  formulation  given  was 


(l+xjy11  +  y1   +\   (l+x)y  =  0     0  .  x  <  1 


ym\   =  y/-n  =  0' 


(IV-9) 


(o)  -  '(1) 


Various  numerical  methods  were  applied  and  close  bounds  for 
the  smallest  eigenvalue  were  given  (  2  ,  p.  236) . 

3.218211  <  \x   <  3.218532 

In  order  to  compare  this  solution  with  the  solution  obtained 
by  the  present  method,  a  solution  was  attempted  by  another  method. 
The  various  steps  involved  were:  1)  The  first  m  terms  and  the 
eigenvalue  were  assumed.  2)  About  40  more  terms  in  the  series  ex- 
pansion were  found  by  using  the  recurrence  relation  obtained  from 
the  differential  equation.  3)  The  first  m  terms  were  then  read- 
justed to  satisfy  the  boundary  conditions.  4)  An  improved  eigen- 
value was  calculated  from  the  Raleigh's  quotient  of  the  polynomial 
solution.   The  cycle  consisting  of  steps  two  to  four  was  repeated 
a  number  of  times,  everytime  improving  the  solution  function  and 
the  eigenvalue. 

However,  it  was  observed  that  the  convergence  was  very  slow 
in  both  of  these  processes.   The  coefficients  in  the  solution 
polynomial  did  not  decrease  rapidly  enough.   And  even  after  about 
20  iterations  the  eigenvalue  agreed  with  the  correct  one,  only  to 
the  second  significant  figure.   To  accelerate  the  convergence,  a 
variable  change  was  made- 
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t  =  2  x 
so  that  the  new  formulation  was 

fcfl+l/2t)£~g  +  2^  +  *(l+l/2t)y  =  0 
dt  a^ 

0  L   t  L   2 

y(o  ■  y(2)  -  °. 

With  this  formulation  agreement  up  to  five  places  was  ob- 
tained in  19  iterations.   More  iterations  deteriorated  the  solu- 
tion rather  than  improving.   The  reason  could  be  attributed  to 
the  fact  that  only  eight  place  arithemetic  was  used,  and  this 
might  perhaps  be  the  closest,  the  answer  can  reach  the  true  solu- 
tion with  the  number  of  places  used  in  the  arithemetic.   The 
solutions  are  compared  in  table  VI. 
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Table  VI 

Results: 

Longitudinal 

Vibrations  of 

a  Cantilever 

Coef f . 
Number 

Comparison 
Solution 

6  th 

Iterations 
13  th 

19th 

X 

3.2184796 

3.2427394 

3.2245237 

3.2185929 

1 

1.0000000 

1.0000000 

1.0000000 

1.0000000 

2 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

■ 

3 

-0.3966840 

-0.3826942 

-0.3909474 

-0.3940352 

4 

0.0661140 

0.0598456 

0.0630494 

0.0642405 

5 

0.0018139 

0.0023530 

0.0020571 

0.0019467 

6 

0.0045958 

0.0040887 

0.0043546 

0.0044540 

7 

-0.0028505 

-0.0024963 

-0.0026736 

-0.0027404 

8 

0.0011162 

0.0009618 

0.0010388 

0.0010680 

9 

-0.0004800 

-0.0004102 

-0.0004452 

-0.0004585 

10 

0.0002170 

0.0001836 

0.0002002 

0.0002064 

11 

-0.0000983 

-0.0000824 

-0.0000903 

-0.0000933 

12 

0.0000445 

0.0000373 

0.0000410 

0.0000425 

13 

-0.0000206 

-0.0000170 

-0.0000188 

-0.0000195 

.  14 

0.0000095 

0.0000078 

0.0000087 

0.0000090 

15 

-0.0000044 

-0.0000033 

-0.0000040 

-0.0000042 

16 

0.0000021 

0.0000019 

0.0000019 

17 

-0.0000010 

-0.0000009 

-0.0000009 

18 

0.0000005 

0.0000004 

0.0000004 

19 

-0.0000002 

-0.0000002 

-0.0000002 

20 

0.0000001 

0.0000001 

0.0000001 

• 

• 

3^ 


iv)  Intermediate  Eigenvalues:-  The  basic  iteration  procedure  al- 
ways converged  to  the  smallest  eigenvalue  and  the  corresponding 
eigenf unction.   However,  the  method  could  be  modified  to  provide 
other  modes  by  orthogonalising  the  initial  trial  with  respect  to 
the  known  modes  and  continually  purifying  the  iterates. 

Thus  the  Iterate  was  required  to  satisfy  all  the  boundary 
conditions  and  the  orthogonality  conditions,  a  total  of  (m+p) 
conditions,  where  the  (p+l)st  mode  was  sought.   This  necessiated 
(m+p)  undetermined  terms  in  the  expansion  polynomial.   Out  of 
these,  m  terms  were  supplied  by  the  constants  of  integration,  and 
the  next  p  terms  in  the  series  were  used  to  provide  the  rest. 

The  orthogonality  condition,  in  the  integral  equation  form, 
cannot  be  directly  used.   The  conversion  to  algebraic  equation 
form  is  illustrated  here.   Even  though  y,  and  y2;  the  first  and 
the  second  eigenf unctions,  are  used  in  the  illustration,  the  pro- 
cess is  perfectly  general.   The  condition  for  y,  and  y2  may  be 
written  as 

\  y2M  [y1l  dx  =  0 


f1 

or   \  y2N[y1]dx  =  0.  (IV-10) 

Both  the  forms  are  equivalent  and  simplicity  will  govern  the  choice, 
For  the  sake  of  illustration  the  later  form  is  chosen. 
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The  first  eigenfunction,  yv ,  would  be  known  In  polynomial 
form.   The  application  of  the  operator  N  will  yield  another  poly- 
nomial 

N[yi]  =  p^x1"1. 

y«  is  sought  as  a  power  series 
y2  =  [ex1"1 

Substituting  in   (IV-10) 


y2N  JyJ  dx  = 


|    (Ic1xJ"l)(Ib1xl-l)dx 
>0   J  l 


■  I   cj  i   Kv1+J"2)to 

j  •'n  i 


i-^i+3-1 


=  LVj  =  o 


where  d3  =  ^I*J=I 
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Thus,  the  conditions  the  iterate  should  satisfy,  are 


BAyl   ■  0    1=1,2 m 


c.dj  =0    1=1,2, 
J 


The  method  was  applied  to  find  the  elgenfunction  correspond- 
ing to  the  second  lowest  eigenvalue  of  the  formulation  (IV-1). 
The  result,  with  the  expansion  of  true  solution,  is  tabulated  in 
Table  VII. 
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■ 

Table  VII 

Results : 

i  Beam  Vibrations  -  Second 

Mode 

..  . 

Coeff . 
Number 

Comparison 
Solution 

Tenth 

Iterations 
Fifth 

Second 

1/A2 

0.00040057 

0.00040205 

0.000509232 

3 

1.0000000 

1.0000000 

1.0000000 

1.0000000 

4 

-2.3561851 

-2.3561802 

-2.3550217 

-2.2766555 

7 

6.9345492 

6.9345261 

6.9090413 

5.4548292 

8 

-7.0024633 

-7.0024360 

-6.9620717 

-4.8558406 

11 

3.4348552 

3.4348092 

3.3694359 

1.1291510 

12 

-2.2072238 

-2.2071817 

-2.1463823 

-0.4514840 

15 

0.3569306 

0.3509095 

0.3287323 

■ 

16 

-0.1681989 

-0.1681846 

-0.1496814 

19 

0.0121331 

0.0121300 

0.0087900 

20 

-0.0045139 

-0.0045122 

-0.0028828 

23 

0.0001725 

0.0001723 

0.0000522 

24 

-0.0000530 

-0.0000529 

-0.0000100 

27 

0.0000012 

0.0000012 

28 

-0.0000003 

-0.0000003 

All  (4n+l)st  and 

(4n+2)nd  coe 

ff.  are  zeros. 

Comparison  solution  = 

cos  3  2 
(cos62x-cosh82x)  +  ST5E3 

-coshS- 

2-sinS2  <sins2x" 

■sinhS2x) 

with  S2  =  \*   =  7.0685465 

obtained  from  10th  iteration. 

■ 

• 

CHAPTER  V 

THE  PINAL  PROBLEM 

Most  of  the  problems  considered  up  to  this  time  were  of 
academic  nature.   It  was  now  decided  to  work  a  problem  of  practi- 
cal interest.   The  problem  selected  is  associated  with  the  name 
of  Graetz[7],  who  obtained  the  first  analytical  solution  of  the 
problem.   The  present  non-dimensional  formulation  was  taken  from 
a  Master's  Thesis  by  Robert  Lipkis  [10]  . 

y11  ♦  yVx  +  (1-x2)  y  =  0     0<.x<_l 

y(0)  =  *(i)  -  °  cr-i) 

These  equations  represent  the  problem  of  'Heat  Transfer  to 
an  Incompressible  Fluid  in  Laminar  Motion' .   A  brief  description 
of  the  problem  is  given  in  Appendix  C. 

Mr.  Lipkis  had  evaluated  the  first  few  eigenvalues  to  a  very 
good  accuracy.   The  present  method  was  applied  only  to  obtain  the 
first  mode.   The  eigenvalue  obtained  (7.31358)  tallys  with  that 
given  in  the  reference  (7.31358)  to  the  given  number  of  places. 

An  interesting  thing  about  this  problem  is  that,  the  quotient 

\      M  y  dx  /  \  N  y  dx 
•'o      '      ^0 

for  an  admissible  function  y  is  also  a  good  approximation  to  the 
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eigenvalue.   In  fact,  when  this  quotient  was  used  by  mistake,  the 
iteration  converged  a  little  earlier  than  it  did  while  using 
Raleigh's  quotient.   This  solution  is  listed  in  column  3  of  the 
Table  VIII. 

The  reference  [1?]  also  listed  the  values  of  y1  at  x  =  1. 
The  value  of  y^,  calculated  from  the  polynomial  solution  by  the 
present'  method  compared  very  well  with  the  value  in  reference. 
However,  this  does  not  mean  much.   Because,  the  mode  function  is 
only  determined  to  within  a  constant  multiple  which  can  always  be 
adjusted  for  any  desired  value  of  yL.   It  only  means  that  both 
solutions  were  normalized  in  some  equivalent  manner. 

A  solution  was  also  attempted  without  the  use  of  Raleigh's 
quotient.   The  term  factored  out  during  normalization  was  taken 
to  be  the  reciprocal  of  an  approximation  to  eigenvalue.   The  pro- 
cess diverged,  oscillating  between  two  values  which  bracketed  the 
true  value.  As  a  remedy  a  weighted  average  of  two  successive 
iterates  was  tried.   This  effected  the  convergence.   But  the  con- 
vergence was  not  fast  enough.   It  is  quite  likely  that  a  proper 
choice  of  weighting  factors  would  do  better.   (The  weighting  pro- 
cedure used  was  nth  approximation  =  0.7  times  the  (n-l)st  approxi- 
mation +0.3  times  the  nth  iterate.)   This  solution  appears  in 
the  second  column  of  Table  VIII,  and  the  solution  by  iteration 
using  Raleigh's  quotient  appears  in  column  4. 
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Table  VIII 

Results:   Heat  Transfer  to  an  Incompressible  Fluid  in  Laminar 
Motion 


Coeff . 

Number 

Comparison 
Solution  1 

Comparison 
Solution  2 

Tenth 
Iteration 

1/AX 

0.13677533 

0.13673189 

0.13073178 

1 

1.0000000 

1.0000000 

1.0000000 

3 

-I.836673O 

-1.8283964 

-1.8283971 

5 

1.3018585 

1.2928571 

1.2928585 

7 

-0.6285^97 

-0.6340986 

-0.6340993 

9 

0.2080135 

0.2202023 

0.2202023 

11 

-0.0536991 

-0.0624800 

-0.0624797 

13 

0.0105566 

0.0143570 

0.0143570 

15 

-0.0017016 

-0.0028671 

-0.0028671 

17 

0.0002170 

0.0004921 

0.0004921 

19 

-0.0000236 

-0.0000758 

-0.0000758 

21 

0.0000020 

0.0000104 

0.0000104 

23 

-0.0000002 

-0.0000013 

-0.0000013 

25 

0.0000000 

0.0000001 

0.0000001 

CHAPTER  VI 
CONCLUSION  AND  DISCUSSION 

Even  though  a  rigourous  proof  of  the  convergence  of  the 
Iteration  for  a  general  problem  is  not  given,  it  could  be  seen 
that  the  process  converges  in  many  cases.   The  method  is  suffi- 
ciently general  in  application,  and  could  be  used  with  advantage 
when  no  other  simple  analytical  solution  is  available. 

The  method  is  equally  applicable  to  equations  with  constant 
and  non-constant  coefficients.   The  only  requirement  is  that  the 
coefficients  should  possess  a  power  series  expansion.   Obviously, 
only  a  finite  number  of  terms  (usually  a  very  few)  could  be  kept, 
so  that  the  error  involved  will  be  dependent  on  the  number  of 
terms  used  and  therefore  on  the  rate  of  convergence  of  the  series 
expansions  of  the  coefficients. 

To  be  able  to  continue  the  iteration  without  getting  stuck 
in  the  integration  process,  no  negative  or  fractional  powers  of 
the  independent  variable  should  appear  in  the  result  of  G[y] . 
The  appearance  of  negative  powers  is  usually  due  to  some  coeffi- 
cient in  the  differential  equation,  when  solved  for  the  highest 
derivative  term.   This  difficulty  can,  in  most  cases,  be  over- 
come by  a  suitable  change  of  variable  e.g.  a  +  t  =  x,  |  a|  >  1. 

Another  difficulty  is  encountered  when  some  coefficient  has 
(a  +  x)  in  the  denominator  with,  |a|  less  than,  equal  to  or 
slightly  larger  than  one.   Division  by  (a  +  x)  is  equivalent  to 
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multiplication  by  its  reciprocal.   The  coefficients  in  the  expan- 
sion of  l/(a+x),  decrease,  remain  the  same  viz.  unity,  or  in- 
crease according  as  |  a|  is  greater  than,  equal  to,  or  less  than 
unity.   It  is  of  utmost  importance  to  see  that  these  coefficients 
decrease  sufficiently  rapidly,  by  effecting  a  change  of  variable, 
if  necessary.   Otherwise,  the  chances  of  convergence  go  down. 

This  case  was  encountered  in  the  problem  of  vibrations  of  a 
tapered  beam  with  a  =  2,  and  in  the  longitudinal  vibrations  of  a 
cantilever.   In  the  .first  case,  the  difficulty  was  solved  by 
switching  the  ends  of  the  beam  and  thus  working  with  a  =  1/2. 
This  in  effect  is  equivalent  to  a  change  of  variable  of  1  -  t  =  x. 
In  the  latter  problem  a  variable  change  of  t  =  2  x,  was  made. 
Both  of  these  changes  sufficiently  accelerated  the  convergence. 

In  the  end,  a  few  words  about  the  use  of  various  subroutines 
would  not  be  out  of  place.   Except  for  subroutines  INTGRD  and 
INITL,  all  subroutines  are  very  general  in  nature,  and  can  be 
used  for  almost  any  problem.   However,  when  the  order  of  the 
equation  is  small  e.g.  2,  the  formation  of  the  D  matrix  is  a  very 
simple  matter.   If  it  is  fed  to  the  program  as  data,  considerable 
computer  core- space  will  be  saved.   This  can  be  used,  with  advan- 
tage, for  carrying  out  more  iterations,  or  for  employing  more 
places  of  arithmetic.   The  same  comment  applies  to  the  use  of  the 
subroutine  BOUND.   For  instance,  the  problem  of  longitudinal  vi- 
brations of  a  cantilever  had  the  following  boundary  conditions- 


(0)  "  *<1) 
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The  condition  yfQ\    =0  leads  to  c„  =  0;  and  then  y/-,\    =0,  means 
simply 

x   i=3 

If  this  is  recognized  in  advance,  a  few  instructions  can  re- 
place the  subroutines  BOUND  and  MATMUL  (A,N,B) ,  and  therefore 
subroutines  DMTRX  and  INVRS  (A,N). 
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APPENDIX  A 

Computer  Programs 

The  programs  for  the  problems  illustrated  in  this  work  are 
listed  here.   Some  of  the  subroutines,  which  are  of  rather  gen- 
eral character  are  listed  in  Appendix  B.  And  Appendix  C  gives 
brief  description  of  the  Final  Problem,  and  the  computer  program 
in  its  entirety. 
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C 
C 


SUBROUTINE    INITL 

THIS    SUBROUTINE     INITIALIZES    THf:    P^CCLC^    OF 

A    BclAK    ON    ELASTIC    FCUNDATION 

COMMON    M  tN i IT i NC iNE  tNP , f PI t NPl f EI 

COKNON    Bt8,8),D(8,U)tC(100)  ,E(10),PI1001 

C(  1)  =  1. 

NC  =  1 

C/'LL    Gf'.TRX 

RETL«iN 

Ei\C 
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SUHKOUTINE  INTGRC 
C      TH  IS  SUBROUTINE  FCRfrS 
C      THE  [NTICRAND  FLR  THE  PROBLEM  OF 
C      A  REAM  CN  ELASTIC  FCUND4TION 

COMMON  MtNf  [TiNCfNEtKPf  KPLiNPl-fCI 

COMMON  P(R ,8) , D(6tl6] ,C(100)  ,6(10)  ,P(  100) 

1  FQRKAT(4E16.81 

2  FQRK/ITC/'i 

WRITEOi  11     fC(  n,I*l»NC) 
KRlTE\3i2J 
C(  1  )=L.-C(  1) 
00    20    I  =  2,:"JC 
20   C(  I  )=-Ct  I) 

WRITEUil)     (C.l  l!tl*ifNC) 

RETURN 

END 
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SUBROUTINE  IN  Ml 
C      THIS  SUBROUTINE  INITIALIZES  THE  PROBLEM  CF 
C      BENCING  OF  A  STkLT 

COMMON  M,N,IT,NC,NE»NP,MPI,NP1,EI 

COMMON  B(P,G) ,0(8, 16) ,C( 100) ,EI10)  ,P(  100) 

E(  1)  =  2. 

E(j2)«~4« 

E(3)=4. 

NE*3 

CALL  DMTRX 

C(  1)  =  1. 

NC«1 

RETURN 

END 
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SUBROUTINE  INTGKC 
C      THIS  SUBROUTINE  FCR>S 
C      THE  INTIGRAND  FOR  THE  PROBLEM  CF 
C     BENDING  OF  A  STRUT 

COMMON  MtN»IT,NC,NE»NP,PPl,NPl,EI 

COMMON    H(%8)  ,0(8,16)  iCCIOOJ  ,  £(  10)  ,  PI  100) 

1  FORMAT  UElV.BJ 

2  FORMAT  (7) 

WRITE!?, 11     (C(  1  ),I  =  1,,\C) 
HRI  r£(3,2) 

CALL    SCiVUL     (CfNCiEfNEfPfNP) 
WRITE (3, 2) 
CI  l)»-4.-«.»P(  I) 
DC    20     1  =  2, TiP 
20    C(  t  )=-4.*P{  I  ) 
NC  =  .\P 

WHITE (3,1)     (C(  1 ),  1*1, NC) 
RETURN 
ENO 
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subroutine  ini  h 
c  this  subroutine   imualizes  the  problem  of 

C      VIBRATIONS  OP  A  BEAP  FIXED  AT  ONE  END 
C      AND  HINGED  AT  ANOTHER  -  FIRST  MODE 

COMHQN  M,Ni  IT,NC»NE|NPiNPl  ,  NP1  ,nl 

CCflKQN  .'}  (  8  ,  8  )  ,  0  (  3  ,  L  6  )  •  C  ( 1 00  )  ,  E  (  10  )  ,  P  (  100  ) 

C(  1)  =  1. 

NC  =  1 

CALL  DfcTax 
RETURN 
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SUBROUTINE  INTGRD 

C      THIS  SUBROUTINE  NORMALIZES  THC  ITERATE  FOR'  THE  PROBLEM  - 

C      VIBRATIONS  OF  A  BEAf  FIXED  AT  ONE  END 

C      AND  HINGCD  AT  ANOTHER  -  ANY  WCDE 
COMMON  MiNi  IT|NCfiNEfNPffPl|NPl  ,EI 
COMMON  P(P,P) ,0(0,16) ,0(100) ,E(  10)  ,P(  100) 

1  FCRMT14E16.8) 

2  FORKATt/1 

DO    20    I«1»NC 

IF(C(  I  )  .•1E.0.0  )    GC    TC    4C 
20    CCM  INUE 
40    EI=C( I  ) 

WRITE (3,1)     EI 

WRITC(3,2) 

DG    60    J«IfNC 
60    C( J)=C( JI/EI 

WRITE (3, 1)     (C(  I  )  ,  i  =  l,NCJ 

RETLRiN 

END 
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subkoutins:  initl 

C      TH[^  SUBROUTINE  INITIALIZGS  THE  PROBLEM  OF 
C     VIBRATIONS  UF  A  rAPEREO  BEAM  of  rectangular 
C     CROSS-SECTIONi  j^'.PLY  SLPPORTEU 

COMMON  M , N , I T , NC »NE  t NP  , PP1 t  NP1 »E1 

COMMON  P.  (  8  ,  R  )  ,  D  (  9  , 1  6  )  ,  C  (  LOO  )  ,  E  (  10  )  ,  Pi  100  ) 
.  2 '  FORMAT ( 14F5.1) 

CALL  DMTRX 

RE4CU,2)  4LPH 

£(  L)-l. 

E(21*2.»I&LPH-1.  ) 

El  5)  =  ( ALPH-L. ) MALPH-l.  ) 

NE-3 

C(  U»l. 

RETURN 

ENO 
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SUBROUTINE  INTGRO 
C     THlj  SUBROUTING  NORMALIZES  THE  ITERATE i 
C      ANO  FORMS  THE  INTIGRANO  FOR  THE  PROBLEM  OF 
C      VIBRATIONS  OF  A  TAPERED  BEAM  QF  RECTANGULAR 
C     CROSS-SECT IOfii  SIMPLY  SUPPORTED 

COMMON  Mt N,  I  Ti  NC  i  NE  , NP  » ivP  L  t  NP  L  •  E I 

COMMON  B( n »  « ) »  D ( ? , 16 ) , C  (  100 )  , E ( 10 )  , P I  100 ) 

1  FORMAT! 4E16. 8) 

2  FORMAT!/) 
NC»40 

C     NORMALIZATION  CF  THE  ITERATE 
00  10  I  =  L,NC 
IF(C( I ) .NE.O.n)  GC  TG  20 

10  CONTINUE 
20  EI-»C(  I  ) 

DO  40  J=IiNC 
40  C(  J  >»C!  J)/EI  .. 

WRITE!3tl)EI 

WRITE! 3i 2) 

WRITE!3ill     (  C  (  i  >  ,  I  =  1  (NC) 
C     '         CONSTRUCTION    OF    THE    INTEGRAND 

NP  =  inC 

CALL    SED IVICNC »E|NEi P , NP) 

WR ITE( 3 i 2) 

DC    50    l»l,NC 
50   C(  I  )*P(  I  ) 

WRITECJ.l)     (C(  I  )  ,1-liNC) 

RETURN 


:iH\j 
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SUBrtOUTfNr  INI  TL 
C      THIS  SUBROUTINE  INITIALIZES  THE  PROBLEM  Or 

C      VI2RATIONS  CF  A  CANTILFVER  BEAM  Or  VARYING  CRCSS-StCT  ION 
COMMON  MtN»ITfNCiNEfNPf^Pl»NPliEI 
COMMON  D(fltS) ,O(0, It) ,C( 100)  , C(10)  ,P (  IOC) 

E(l)*2. 
Ef 2)»l. 

NE  =  2 

CALL  DMTRX 

C  (  1  )  =-  '» . 

C(2)=0. 

C(3)=l. 

NC  =  3 

RE?  TURN 

END 
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SUBROUTINE  DMTRX 

c 

THIS  SUBROUTINE  CLNSTRUCTS  D  MATRIX  FROM  E  MATRIX 

c 

FOR  THf-  PROCLEf!  CF 

c 

VIQR4T10NS  CP  A  CANTILEVER  BEAM  OF  VARYING  CROSS-SECTION 

DIMENSION  F(10) 

COMMON  K,N|lT|NC|NE»NPtfPl»NPl,El 

COMMON  B  (  ?, ,  C  )  ,  0  (  8  ,  It )  ,  C  (  100  )  ,  E  (  1  0  )  ,  P  {  100  ) 

1  FORMAT (4E16. 8) 

2  FORMAT(  14F-5.2) 

F(  1  )=l. 

DC  20  1  =  2,  M 

20  F(  I  )  =  F (  I-1)«FLGAT(  I-L) 

REAC(li2)  (  16(  I  >J)i  J»l  iM  i'1-liMJ 

DO  100  1=1, H 

IF(  I.GT.N)  GO  TO  iO 

OU  AG  JaLtM 

4(J  D(  I  ,J)=B( I ,J)*r (J) 

GC  TO  ICO 

• 

60  D(  I  ,1)=S (  I i i) 

DC  80  J  =  2,M 

• 

D(  I  ,J)=C(  1,1) 

, 

DC  £0  K  =  2,J 

jk  =  j-:; 

80  D(  I  ,J)=0(  I  ,  J)  +  F(  J)*f»  (I  |K)/F(  JK  +  1) 

100  CONTINUE 

DC  120  I  =  1 ,  M 

C0N*1. 

3 

00  120  J«1#M 

0  (  I  ,  J  )  =  D  (  I  i  J)»C'ON 

120  C0fci?CQN*2. 

CALL  IflVRSlOfM) 

• 

MR  I  TE  (3. 1 )  (  {  0  (  I  ,  J  )  ,  J  =  1 ,  M  )  ,  I  =  1 ,  M  ) 

RCTURN 

E.NG 

» 
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SUBROUTINE  IN1GRC 

c 

THIS  SUPiUJUTlNr:  NCRPALUES  THE  ITERATE, 

c 

FINCS  ITS  RALEIGHS  CUCTIENT  ANO  FORMS 

c 

THE  [NTIGPANO  FCR  THE  PROBLEM  UF 

c 

VIBRATIONS  UP  A  CANTILEVER  BEAM  OF  VARYING 

DIMENSION  F( 1001 p G 1  ICC )  ,R( 100) 

COMMON  M,N»IT,NC»hE,NP,FPl,NPi,EI 

COMKON  I*  (8,8) ,  D(8, 16)  ,C  ( 100)  ,E  (10)  ,P  ( 100  ) 

CROSS-SECTION 

1 

FORf  AKAE16.8) 

4 

fckmti/  ) 

c 

NORMALIZATION  OF  THE  ITERATE 

DC  LO  I=1,NC 

IF  (  C(  I  )  .NCO.O)  GC  TO  20 

IC 

CGNTINUF. 

• 

20 

EI*C( 1 ) 

DC  <iO  J=I  iNC 

40 

C( J)=C( J) /EI 

WRITE!  3,  DEI 

WSI TE(3,4) 

WRITE! 3, 11  (C( I )  ,1  =  1 »NC) 

c 

t 

CALCULATION  OF  RALEIGHS  QUOTIENT 

CALL  DIFFER  ( C  ,  NC  ,  F ,NF  ,  1  ) 

WRITE( 3, A) 

CALL  OIFFFR  (F  ,!\F  ,G,NG,  I) 

WRITE (3,4) 

CALL  SEMUL  (G»NG,£,NE»R,NR) 

W«  ITE(3,4) 

DO  bO  I=1,NF 

50 

R(  I  )=P( I  )+F(  I ) /2. 

CALL  S-EMUL  (C,\C,k,.\R,G,NG) 

WRITE (3, A) 

CALL  INTGRA  (G,NG,R,.\:"<  •  1) 

Wrt I TE( 3, A) 

DM  =0. 

CG\'=1. 

DO  60  I«3,NR 

, 
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CCN=CCN*2. 
60  DM*CM4C0N*ft( 1-1 ) 

CALL  SEHUL  (CiNCiEfKEfG'iNG) 

WRITE  13, 41 

CALL  SEKUL  (  C  ,  NC  t  G  ,  ,\G  ,  R  ,  NR  ) 

WRITE! 3,4) 

CALL  INTGRA  IR,NR»G,NG,  I) 

WRITE (3, 4) 

DN  =  C. 

CCN«1. 

DC  70  I=3,NG 

C0N=CUi\*2. 
70  DN«CN-CO.N»G(  I-  1 ) 

EI*4.»DM/0N 

WR  HE[3,1)  EI  t  OK  t  UN 

WRITE (3,4 ) 

CONSTRUCTION  OF  THE  INTEGRAND 

CALL  SCDIV  IF,NF,E,NE,G,NC) 

WRI  IE  (3,4) 

DO  60  1  =  1,  'iC 
80  C(  I  )=-FI*C(  I  )/4.-G(  N./2. 

WRI  TO (3,1)  (C(I)tl-lfNC) 

WRITE (3, 4) 

RE  TORN 

ENO 


59 


SUP ROUT INC  BOUi-iO 
C      THIS  SURROUTINfc  CCMPUTES  THE  F  VECTOR  FCR  THE  PROBLEM  OF 
c     vidrWioms  OF  A  CANTILEVER  BEAM  OF  VARYING  CRCSS-SECTION 

DIMENSION  Y( 1C  ) 

COMKCN  M,Ni  [T|NC»NE»NPfPPiiNPl,EI 

COKtfCN    R(8,8)tO(8,U)fCllOO)  ,E( 10) ,P< 100) 

1  FORMAT ( AC  16.8) 

2  FORMAT!/ ) 
00     10     1=1, M 

10   C(  I  )=C. 

C0N*2*»»H 

DO    30    [*iiW 

Y(  I  )=0. 

A  I  =  I 

DO    20    J*MP1»NC 

Y(  I  )  =  Y(  I  )+CCN*C( J) 
2  0    CCN=CO,\«(  FLOAT  I  J  )  /  (  FLO  A  T(  J  )  +  1 . -AI  )  ) '2  . 

C0N=2.*«N 
DO    20    J=l,I 
30    CON«CON»  FLOAT  (  ,\+l-J)/2. 
HRITE'('iil)     CC(I'»tl"l.M) 
RETURN 
ENO 
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C 

c 
c 
c 


MAIN  PROGRAM  FORR  OBTAINING  THE  SECOND  • 

EIGCNFUNCTION  FOR  THE  PROBLEM  Of: 

VIBRATIONS  GF  A  BCAf  FIX6D  AT  ONE  END 

AND  HINGED  AT  ANOTHER 

DIMENSION  RUOO) 

COMMON  M ,N , IT , NC , NE ,NP i fPl »NP1 ,EI 

COMMON  Rt«f6),D(fi',  U.)fC(lOO)  ,5(10)  ,P(100) 

1  FORMAT (4E16. 81 

2  FORMAT  1 14151 

3  FORMAT  I  IHi'l 
A  FORMAT! /) 

R E A C ( 1 1  2 )  M t N i "IT L 

MP1=M+1 

NPl*N+l 

CALL  INI T  L 

WRITE(3»3] 

DC  40  IT=1»ITL 

CALl  INTGRO 

CALL  I NTGR A(C i NC  fR i NR , M ) 

DC  20  I=MPl,NC 
20  C(  I  1*R(  I  ) 

CALL  BOUND 

WRITE I  3,4] 

CALL  MATMl)L(DiMPliC) 
AC  WRITE (3, 3) 

CALL  RESULT 

STLP 

ENC 
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SUDrtOUTINC  INI  1 L 
C      THIS  SUBROUTINE  INITIALIZES  THE  PROBLEM  CF 
C      VIBRATIONS  OP  A  BEA/'  FIXED  AT  ONE  END 
C      ANC  HINGED  AT  ANOTHER  -  SECOND  KOOE 

DIMENSION  RdOGl 

COMMON  M,M,  ITiNCiNE  »NPf  PPliNPl  ,EI 

COUPON  C  (  8  ,  e  )  ,  D  (  B',  1  6  )  ,  C  (  100  )  ,  E  (  10  )  ,  P  {  100  ) 

1  F0RPAT<4E16.8) 

2  FORMAT { 141 5 J 
REALI1.2)    NP 
REACH, 11 CR< I )  »  1=1, NP) 
DO    20    I=1,NP' 

A  I  =  1 
P(I)=U. 
DO    2C    J=1,NP 
20    P(  i)aPtI)+R(JJ/(AI*FLOATIJ)-l.l 
WRITE13,1J     (P(  I  )  ,1=1, NP) 
CALL    DMTRX 
C(8)*l« 

nc=£ 

CALL    BOUtJU 

CALL   MATMUL(D,MP1 ,C) 

CC7)»C.(5) 

C(  fiJ  =  0. 

RETURN 

ENO 
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C 
C 
C 


SUBROUTINE  DMTKX 

THIS  SUBHCUTINt:  CONSTRUCTS  D  MATRIX  FCK  Trill    PROBLEM  OF 

VIBRATIONS  OF  A  BEAM  FIXED  AT  ONE  END 

AND  HINGED  AT  ANCIHCR  -  SECOND  MODE 

DIMENSION  F(IQ) 

COMMON  MiNtITfNCtNEfNP,MPlfNPl,EI 

COMMON  B(8f8)»D(8f  l6)fC(100)  fE(l.O)  tP(lOO) 

1  FORMAT (4E16. 8) 

2  FQRMATt 14F5.2) 

3  FORMATI/1 
F(  1  )  =  1. 

DO  ?.C.    I»1,MP1 
20  F( I )=F( I-l)»FLUA Hl-l) 

REAC(lt2)  I(B( I,J),J*l,M),I*liM) 

DO  ICO  1  =  1, M 

IF(  l.GT.M)  GO  TO  6  0 

DO  40  J»1,MP1 
40  D(  I  ,J)=R(  I  ,J)*F*U) 

GO  TO  100 
6G  DC  I,1)=BU,1) 

DC  60  J»2,MP1 

D(  1  ,J)=r (  I  ,  I) 

DC  SO  K«2,J 

JK^J-r; 
6C  D(  I,J)=D(I  ,  J)+FU)*P.  (I  ,Kl/F{ JK+1) 
ICO  CONTINUE 

DO  12C  1  =  1, MP  I 

120  d(mpii  i )*?n ) 

WRITE! 3,1)  {  (D(  I  ,J)  ,J*l,MPll  f  1=1  f  MPD 

CALL  INVRS(D,MP1) 

WRITE! 3,3) 

HR I TE ( 3, I J  ( ( C (  I  ,  J  )  , J? 1 , MP  1 )  , I  =  1 »  M P 1 ) 

RETURN 

ENO 
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SUBROUTINE  BOUNC 
C     THIS  SUBROUTINE  COMPUTES  THE  F  VECTCK 
C      IN  1H1-  MATRIX  ECUAHGN  CC»F  FOR  THE  PROBLEM  OF 
C      VtDRATIUNS  UF  A  PLAi*  FIXED  AT  ONE  END 
C     AND  HINGED  AT  ANOTHER  -  SECOND  MODE 
D  I  KENS  ION  Y(1G) 

COMMON  M,N»IT,NCiNCfNP»*PliNPitEI 
COMMON  B ( 8 , F ) , C ( R  ,  1 6 ) , C ( 100 ) , E 1 10 )  , P (  100 ) 
1  FCRMAn'iE16.8) 
DO  20  I=1,MP1 
20  C(  I  )=C. 
CCN=1. 
DO  60  1=1, M 
Y(  I  )  =  0. 
AI*I 

DO  40  J=8,NC 
Y(  I  )=Y(  I )+CON»C(J) 
AG  CON-CON* -(FLOAT! J ) / ( FLOAT ( J )  +  1 . -A  I  )  ) 
CON=l. 
DO  60  J=l, I 
60  C0N=CQN*(8.-FL0ATIJ) J 

DC  80  J=1,M 
8  0  C(  I  )=C( I )-B ( I , J)«Y( J) 
DO  ICO  I=8,NC 
ICO  C(MPl)=C(MPl)-P(  I  )*C( I  ) 

WRITE  (3,1)     (C(  I)  ,  l  =  l,.v?l) 
WRITE!  3i  1.1     (C(  I  )  i  1=1.  NO 
RETURN 
END 
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SU8UQU1 INF  INI J'L 
C      THIS  SUBROUTINE  INITIALIZES  THE  PROBLEM  CF 
C      HCAl  TRANSFER  TC  AN  INCCRPRE SS I  RLE  FLLIC 
C      IK  LAMINAR  MOTICN 

CONM3N  M.NtlTiNCtNetNPt^PltNPlfEI 

COMMON  BC  8  »8> , 0(8,16) ,C ( 100) ,E ( 10) ,P< 100) 

CALL  DMTRX 

G(l)=-1. 

E12)=C. 

E  (  3  )  =  1 . 

NE  =  J 

C(  1)  =  1. 

C( ?)=C. 

C( j)=-1. 

NC»3 

RETURN 

ENO 


. 
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' 

SUBROUTINE  INTGRC 

c 

ThIS  SUBROUTINE  NORMALIZES  THE  ITERATGi 

c 

FINOS  ITS  UALEICHS  QUOTIENT  AND  FORMS 

c 

THE  INTIGRANO  FOR  THE  PROBLEM  UF 

c 

HEAT  TRANSFER  TC  AN  INCCMPRESSIBLE  FLUID 

c 

IN  LAMINAR  MOTION 

DIMENSION  F(lOC)iGllOO)  ,H(100)  ,R!100) 

COMMON  M»N»lT,NC»NEfNP,MPl,NPl,£I 

CO»'CN  Bl8,e),D!8,U),C!l00)  ,EI 10) ,P! 100) 

L 

FORMAT (4E16.0) 

2 

FORMAT!/) 

c 

NORMALIZATION  CF  THE  ITERATE 

DC  20  1*1, NC 

IFICC I J.NE.0.0)  G'C  TC  3C 

2  0 

CONTINUE 

30 

CIJ=C( I ) 

DO  40  J«1,NC 

40 

C ( J)=C( J)/CIJ 

WRITE!  3,1)  (C(  I  )  ,  1=1,  NO 

WRITE! 3, 2) 

c 

CALCULATION  OF  RALEIGHS  QUOTIENT 

CALL  SEMUL I C  »NC , E ,NE ,G , NG ) 

WRITE!3,2) 

CALL  SEMUL (C iNC f GCN6 1 H v NH1 

WRITE! 3, 2) 

F(1)=C. 

CALL  lNTGRA.(HfRH|F>KFf  1) 

WRITE! 3,2 1 

DN  =  C. 

DC  80  I*l,NF 

• 

ec 

DN»CN+F! I) 

CALL  OIFFERtCt NC,F,NF, I ) 

WR  I  fE!3,2) 

CALL  niFFER!F,NF,h,KH, I) 

WRlTE!3i?) 

DU  60  1  =  1  ,NH 
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60    H(  I  )=H( I  )+F( 1  +  1) 

WRITEJ  3i  1)     (IK  I  )  »  I  =  ifNH) 

wi>  I  r R  (  3,  ?) 

CALL    SEMiH.     (C,NCiH,NH,R,NR) 

HH  I  TECi,2) 

H(  1  )=C. 

CALL    INTGRA    |«»iNRtH,NH,l) 

HRI  r  E ( 3 , 2 ) 

D*  =  C. 

L)G    70    I«lvNH 
70    DM»UM+H( I ) 

H(  L)=G. 

E [=CM/0N 

WR'lTCOi  I)    Eli  DK.fCN 

MRIT£(3>2) 

CONSTRUCTION   OF    THE    INTEGRAND 

NC*NG 

no  90  i*ifNC 

C(  i  )=r,(  i  )-F(  i  +  i)/ai 

IF(  I.CT.MF)       C(  1)=G(  I) 

go  CONTINUE 
RIU.»0. 

WRITE  I. 3,  I)     (C(  I  )  ,  1  =  1, NO 

WRI  fC('3,2) 

RETLr<.s 


67 


APPENDIX  B 

Programs  and  Subroutines  of  General  Nature 

The  main  program  and  all  the  subroutines  listed  in  this  ap- 
pendix are  very  general  in  nature,  except  that  the  subroutines 
DMTRX  and  BOUND  can  only  be  used  if  the  boundaries  are  x  =  0  and 
x  •■  1. 
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C     MAth  PfUlGRAM  FOR  ALL  EQUILIBRIUM  PROBLEMS  AND  FOR 
C      THE  FIRST  MOOE  L  F  ALL  EIGENVALUE  PROBLEMS 
C      M  IS  THf  ORDER  CF  CCUATICN  AND  N  IS  THE  NUMBER 
C     OF  BOUNOA«*Y  CONDI TICNS  'U  X=0 
DIMENSION  R( LOO) 

COMMON  MfN|  IT,  N'C  ,i\E  ,  NP  ,  NPl »  NP1  •  E  I 
COMMON  B ( S  ,  8 )  , 0 ( 8  ,  16 ) , C ( 100 )  , E ( 1 0 )  ,  P (  100 ) 

1  FGKrAT('.E16.0) 

2  FQRM4T(SI5) 
.i  FORMAT  J  IH1  ) 
A  FOKMATt/1 

REAC(1V2)  MtNi ITL 

MP1*M+1 

NP  L»N+1 

CALL  INITL 

WRITE(3,3) 

DG  40  I  T  =  1  ,  ITL 

CALL  INTGRC 

CALl  INTGRA(CiNCiKiNR,M) 

NC=NR 

DG  20  I«MPltNC 
20  C(  I  )  =  '•'(  I  ) 

CALL  BOUND 

CALL  MATKUL  ID»M»C) 
AG  WR  I  TE  (  3,  3) 

CALL  RESULT 

STOP 

END 
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SUBROUTINE  DMT ft X 
C      THIS  SUBROUTINE  CCNSTRUCTS  D  MATRIX  F.<C,V  B  MATRIX 
C      FC'i  ALL  EQUILIBRIUM  PROBLEMS  AND  FOR  THE 

c    fikst  hooe  cf  all  eigenvalue  problems 

DIMENSION  FC 10) 

COMMON  M,N»IT,NC,NE,NP,MP1,NP1»£I 

COMMON  Q(B,  fl) ,0(6,16) ,  C (100) , E( 10) ,P( 100) 

1  FORMAT (AE16. 8) 

2  FORMAT! UF5.1) 

3  F0RKATC//20X,  tHG    MATRIX  /) 
F(  1)  =  1. 

.  DG  20  1  =  2, M 
2  0.  F(  I  )  =  F(  I-l)»FLL,vr(  1-1) 

REAL.  (1,2)  (  (B(  I  ,J)  ,J=1,MP1)  ,  I=1,M) 

DC  100  1=1 tM 

IF(  I.GT.N)  CO  TC  .60 

DC  40  J«1,M 
AC  D(  I  ,J)»B(  I  ,J)  »TW) 

GC  TO  ICO 
60  D(  I  ,1)=[J.  (  I  ,1) 

DG  80  J=2,N 

D(  I  ,J)=F.  (  I  ,1) 

DC  CO  K=2,J 

j  k  =  j-;-; 

8  0  D(  I  ,J)=D( I  ,J)  +  F( J)«e (I  iK)/F( JK  +  l) 
100  CONTINUE 

CALL  INVRS(OtM) 

HKIIE(3,3) 

WRITE!  3,1)     (  (D  (  I  ,  j)  ,  J=l  »M)  ,1  =  1  ,M) 

RETURN 

END 


• 
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c 
c 

c 
c 

SUBROUTINE  INVRSCA.N) 

ThIS  SUBRCUTINL  INVERTS  MATRIX  A  OF  ORDER  N 

BY  GAUSS-5  lEDIiL  RGDCdlCN 

DIMENSION  A  18, 16) 
2  FORMAT (2{//20X,lCHNC  INVERSE  //)) 

AUGMENTING  THE  MTRIX  A  PY  AN  IDENTITY  MATRIX 

NN=N+N 

DC  20  1=1, N 

IN*I+N 

DO  10  J=1,N 

JN-J+N 
1U  At  I  ,J,\)  =  0. 
2  0  A  (  I  ,  I N  1  =  1 . 

THC  REDUCTION  PROCESS  STARTS  HERE 

DC  ICC  M>lfN 
30  DIV*A(M»M) 

IF (CIV. CO. 0.0)  GC  TC  7  0 

DO  40  J=1,NN   • 
40  A(M|J)*A(H,J)/DIV 

DC  60  1=1, N 

If II.EQ.M]  GO  TC 6  0 

A  I. '■  =  .'(  I,M) 

DC  tO  J=1,NN 
50  A(  I  ,J)=At  I  ,  J)-AIK«A(M,  J) 
60  CONTINUE 

GC  TO  ICO 
70  DC  SO  I  =  M,N 

- 

IFIAdtMJ.EO.O.O).  GC  TO  <30 

c 

DC  LGOP  80  EFFECTS  INTERCHANGE  OF  ThC  ITH  AND 
■DC  60  J=l, NN 
DUMY*A(  I  ,  J  ) 
.  A(  I  ,J)=A(K,J) 
60  MM,J)=DUKY 

GC  TO  30 
90  CONTINUE 
WRITS! 3.2) 
GC  TO  120 

MTH  ROWS 
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100  CONTINUC 

TRANSFERING  INVERSE  CF  A  IN  THE  PLACE  OF  A 

DC  110  [»l,N 

DC  110  J  =  l,N 

JN*J+N 
110  At  I  ,J)  =  A(  I  ,  J.\) 
120  RETURN 

END 
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SUflKQUTINE  BOUi.C 
C      THE  SUBROUTINE  COMPUTES  THE  F  VCCTUx 
C      FGK  ALL  EQUILIBRIUM  PROBLEMS  AiiO  FOR  THf- 
C      FIRST  .''CDC  OF  ALL  EIGENVALUE  PROBLEMS 

DIMENSION  Y( 10  ) 

COMMON  MtN  i  IT,  NC  i  M:  ,NP  ,  f;P  1 ,  NPl  ,  E  I 

COMMON  B(8,R)  ,0(8,  16)  ,C{1.00)  ,E(  10)  ,P(  100) 
"  1  FORMAT (AC  16. 8) 

2  F0RMAT(//20X,5HB0UNC/J 

DC    10    I  =  1,M 
10    C(  I  )  =  ('.(  I  ,NP1) 

CGN*l. 

DC    30     I«1,M 

Y(  I  )  =  0. 

A  I  =  1 

DC    20    J! M=M Pi,.\C 

AJM*JM 

Y(I)=Y(I  »+CON*'C(Jh) 
20   CON*CGN«(AJM/(AJK*lv-A I  )  ) 

CC.\  =  1. 

DO    30    J=li I 

AJM*M*l-J 
30    CCN=CONVAJM 

DC    AG    I«NP1,M 

DC    40    J=1,M 
4  0    C(  I  )=C(  I  )-p,(  I,  J)«V(J) 

WRITE! 3,21 

WRITE!  3,1  MCI  II,  1*1, Ml 

RETURN 

END 
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SUBROUTINE  KATHUL  IA,N,N 
C      THIS  SUBROUTINE  MULTIPLIES  A  SQUARE  MATRIX  A 
C    ■  OF  CRIJGR  N  BY  A  COLUMN  VECTOR  8.  THE  ANSWER 
C      APPEARS  IN  THE  COLUMN  8. 

DIMENSION  A { 8 , 16 ) , B { 8 ) , 0 ( 8 ) 
1  FORMAT (4E16.8) 
3  FORMAT  (20Xi6H.yAT.MLL) 
WRI  fE(3,3) 
DO  100  1=1, N 
D(  I  )  =  0. 
DG  100  J  =  1,N 
100  D(  I  )*DI  I )+A( I , J) »B(  Jl 

DG  20C  1  =  1,  N 
200  B{  I  )=D(  I  ) 

WRITE (3,  1)  (H(  I  )  ,  1  =  1  ,N1 

RETURN 

END 
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SUBROUTINE  SjEMUL  [A,NA ,H,NB,CtNC) 
SERIES  A  *  SERIES  B  =  SERIES  C 

DIMENSION  A( 100)  ,S ( ICC)  ,Cl  100) 
1  FCUi'.ATI  4E16.8) 

3  FORKATI/ZOX.ZIHSEkIES  multiplication  /) 
WRITE! 3, 3) 
NON&+NB-1 

IF(NA.GT.MR)  GC  TC  101 

DG  IOC  1=1, NC 

MI=I 

C(  I  )»0. 

IF{  I.GT.NA)   MI  =  .\A 

DG  100  J=l ,MJ 

MIJ*I-J+1 
100  C(  I  )  =  C(  I  )+A(  J)  *Mi«  I  J) 

GG  TG  201 
10  1  OG  200  I  =  l,i\C 

M  I  =  i 

C(  I  )  =  C. 

IF(  l.GT.NP  )   MI«NB 

GC  200  J  =  1,M 

MIJ=I-J+l 

200  C(  I  )=C(  I  HSI  J)  »A(MIJ) 

201  WRITE! 3,  I J  (C(  I  )  ,  1  =  1,  NO 
RETURN 

END 
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SUBROUTINE  SGOIV  {  A  ,  NA  ,  I!  ,fi'd  ,  C  ,  NC  ) 

DIMENSION  A(  LOO  ,  B  {  100  )  , C  (  100  ) 
C     SUBROUTINE  F(j:<  DIVISION  OF  TWO  SERIES 
C     SERIES  A/SCRIES  B=SERIES  C 
1  FORMAT  UE 1 6.  81 
3  F0RMAT(//20X,6HSEDIV  /) 

WRITE (3, 3) 

1  =  ] 

C(  I  )=A(  I  )/0(  I) 
10  1=  Ml  ■ 

M I « 1 

IFU.GT.NB1       M1«NQ 

C(  I  )=A(  I  ) 

DC  20  J=2,MI 

MIJ=I+1-J 
20  C(  I  )=C(  I  )-R(  JI*C<*  Ul 

C(  I  J=C(  I  1/8  I  1) 

IF (  I.LT.NC J  GO  TC  1C 

WRITC (3,1)  (C(  I  )  ,  I=1,NCJ 

RETURN 

END 
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SUBROUTINE    DIFFEtUG.NC  ,CtND|N) 

NTH  DERIVATIVE  LF  SERIES  C  ■  SERI.ES  D 

DIMENSION  C(IOO)  id  IGO) 

1  FQRKAT(4E1'6«8) 

3    FORM!     (/2CX,     IC    HDERIVATIVE       /) 

WRITE! 3* 3) 

AF  =  1. 

DO     10    1=2, N 

A  I  =  I 
10    AF  =  AF-*A1 

iFU.EG.l)       AF  =  1. 

DU     ?0     I=1,NC 

IN=1*N 

AI  =  1 

A  I  H  =  I N 

[)(  I  )  =  C(  IN)*AF 
20    AF=AF»AIN/A] 

■NO»NC-N 

WKIT5(3,1)     (0(11* I"    I i NO) 

RETURN 

END 


78 


SURKOUTINE    RESULT 
C  -THIS    SUBROUTINE:    GIVES    THE    RESULTS    Cf:    AN 

C  E0U1LIBR IUM    PROBLEM 

CCMM3N    M,N|  IT,i\C,i\E  ,KP,P?l,HPl  ,El 

COMfON    B(8«e)f  0(8fl6)«.C(l00)  tECIO)  tPClOOJ 
1    FORfcAT(4E16.8) 

HRMEOtll     (CU  It  I- If  NC) 

RETURN 

ENC 
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SUBROUTINE  RESULT 

C      THIS  SUPROUTINf:  GIVES  ThE  RESULTS  OF  AN 

C     EIGENVALUE   PRGPLCM 

COMMON  M,Ni IT|NCiNE,NP|NPliNPl ,EI 

COMMON  R  (Ri8  ) 1 -0.(6, 16  )  ,C(  100)  ,E(  10)  ,P(  1001 

1  FCRPATUE16.8) 

2  FQR?AT(/1 
WRITE! 3  til     FI 
WRIT£(2,1.1    FI 
WRITEOi?) 

OC    20    [«1,NC 

IF(CU).NE.O.O)    GC    TC    4C 
2  0   CONTINUE 
AG    EI=G( I ) 

DC    60    J»l |NC 
6  0    C{ J  )=(..( J) /EI 

WRITE  (3,1)     (G(  I  )  ,  1=1,  NO' 

WRITE12.1)     (C(  I  )  ,  1=1, NO 

RETURX 

END 
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APPENDIX  C 

• 

Description  of  Final  Problem 

This  is  a  steady  state  heat-transfer  probl 

em.   At  steady 

state  there  are  no  variations  with  time  and  the 

energy  equation 

for 

a  fluid  in  laminar  motion  in  a  cylincrecal 

tube  can  be 

written 

as 

(1) 

where  T 

=  temperature; 

a 

=  radial  co-ordinate  measured  from  the 
tube  radius 

tube  axis;  a  , 

U  =  mean  flow  velocity 

m 

o 

=  thermal  diffusivity, 

z 

=  anial  co-ordinate  measured  from  the  e 
pipe. 

ntrance  of  the 

The 

equation  reflects  the  following  postulates: 

1) 

The  fluid  is  incompressible; 

2) 

No  heat  is  transferred  by  conduction  ir 

i  the  z  direction. 

3) 

The  velocity  distribution  is  parabolic 

at  all  cross- 

sections  and  is  maintained  independent 

of  the  temperature 

(i.e.  viscosity  and  density  are  not  the 

functions  of 

■ 

temperature) ; 

4) 

Thermal  conductivity  and  the  product  y0p  are  independent 

of  temperature. 
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The  boundary  conditions  specified  were  that 
1)  the  fluid  and  wall  were  uniformly  at  temperature  T,  for 
z  ,<  0, 

i.e.   T  =  T,',    z  <,  0 

and  2)  the  wall  temperature  varies  linearly  for  z  >  0. 


T  =  T^l  +  -2-z)  a  =  aQ 

o 

The  following  dimensionless  variables  were  defined. 


X    = 

a 
ao 

U)  = 

Z      _a 

a  2   2U 
o          m 

t  = 

T-^ 
Ti 

In  terms  of  these  dimensionless  variables  the  formulation  (1) 
would  be 

X  3X  VX3x'     U  X  ;3ui 

A  product  solution  of  the  form 

t  =  1  -   (x)    U) 
was  employed  to  separate  the  variables.   It  was  found  that 

W  =  e 


82 


2 

where  -g  was  the  constant  of  separation.   The  equation  for  the 


y  was 


&  +  1  |X  t   62(1.x2)y  .  0 

QX 


with  the  boundary  conditions 


€  -  «     at  x  -  o 


y  =  0    at  x  =  1. 


The  computer  program  in  its  entirety  is  given  in  the 
following  pages. 
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C 

c 
c 
c 


MAIh  PROGRAM  FOR  ALL  EQUILIBRIUM  PROBLEMS  *NO  FCR 

THE  r  I R T, T  MODE  LF  ALL  EIGENVALUE  PROBLEMS 

M  IS  fHf.  ORDER  Cf  EQUATION  AND  N  IS  THE  NUMBER 

OF  BOUNDARY  CONDITIONS  AT  X=0 

DIMENSION  R'flOC) 

COMMON  M ,N , IT, NC ,NE ,NP  »  MP1 ,NP1 iEI 

COMMON  RCft.«)tD(8,UJ,C{ipO)iC(10).P(  100) 

1FCRMAT{4E16»8) 

2  FORMAT (5  IS) 

3  FORMAT ( 1H1 ) 

4  FORMAT!/ 1 

R 6  A C 1 1 , 2  )    M t N t IT L 
M?l=M+l 
NP  l.*N*l 
CALL  I  (ill  I. 
WRITE  (3,  3) 
CC  40  IT«1,ITL 
■  CALL  INTGRO 
CALL  INTGRA(C,NCik»MRf M) 

NC=KR 

DC  20  I*MPlfNC 


20  C(  I  >-R( I ) 
CALL  COUND 
CALL  MATMUL 

AG  WP  ITE(3,3) 
CALL  RESULT 
STOP 
END 


{  D  1 1-  i  C ) 
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C 
C 

c 


SUBROUTINE  INITl 

THIS  SUBROUTINE  INITIALIZES  THE  PROBLEM  OF 

HEAT  TRANSFER  TC  AN  iNCCNPRESSIBLE  FLUIC 

I,\  LAMINAR  MOT  ION 

COMMON  H»N»-ITiNCtNEfNPi*'Pl»Npl  'f  I 

COMMON  B < 8 ,8 ) , V i 8 , U ) , C ( 100) , E ( 10 ) , P ( 100 ) 

CALL  DMTRX 

Ell)*-  1. 

E( ?)=0. 

E( 3)=l. 

NC  =  3 

C(  1  )  =  l. 

CI2J*0. 

C( 3)=-l. 

NO  3 

RETURN  • 

ENU 
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C 
C 

c 


SUBROUTINE  DMTKX 

THIS  SUBROUTINE  CONSTRUCTS  D  MATRIX  FROM  B  MATRIX 

FCk  ALL  EQUILIBRIUM  PROBLEMS  AND  FOR  THE 

FIRST  MODE  OF  ALL  EIGENVALUE  PROBLEMS 

DIMENSION  FUO) 

COMK.N  M,N,ITfNC,NE,NP,.*Pt,NPl,EI   . 

CCMKCN  B*B,8) , D(8 , 16) ,C (100) ,E ( 10) ,P<  100 ) 

1  FORMAT (4E16. 8) 

2  FORMAT \ 14F5.1) 

3  FORMAT  (//?.0X,  BHU  MATRIX  /) 

F(  l)»l. 
DO  20  I  =  2,iv. 
2  0  F(  I  )  =  Fl 1-1 )»FLOA T( 1-1) 

REACClt?)(CB(ItJ»tJ-lt»*Pl)iI*ltHJ 
DC  100  I»liM 
IF(  I.GT.N)  CO  TO  60 
DG  AC  J=liM 
AG  D(  I,J)*BUtJ)*F(J) 

GO  TO  ICO 
60  D(  I,l)=BU.l) 
DC  80  J  =  2,M 
Dt  l,J)=B(  I  il  ) 
DC  eo  K=2,J 

JK»J-K 
80  D(I,J)=DU,J)  +  F{J)*B(I,K)/F{JK+1) 

100  CONTINUE 

CALL  INVRS(OtK) 

WR  ITC(3f3) 

WRITE  I  3,1)  I  IDC  It J)i .J«liK)iI*l»K) 

RETURN 
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SUBROUTINE  INVKMAth) 
C      THIS  SUBROUTINE  INVERTS  MATRIX  A  OF  ORDER  N 
C     HY  GAUSS-SIEOEL  KEDUCJICN 

DIMENSION  MB,  lc) 
2  FORMAT! 2 1//20X,1CHNC  INVERSE  //)) 
C      AUGMENTING  THE  ,V A T K 1  X  A  BY  AN  IDENTITY  MATRIX 

NN*N*N 

DO  20  I»l,N 

IM=  I  +.\ 

DO  10  J=1,N 

i&  a(  I  ,jn)=g. 

2  0  A  (  I ,  IN )  *  I . 

c    the:  reduction  process  starts  here 

DO    100    M«1,N 

3  0    DIV«MM,M) 

IF(CIV.EQ.O.O)     GG    TC    7  0 

DC    4  0    J=1,NN 

ao  AC-, ,j)=a{'-' ,j )/;:; V 

DO    60    1=1, N 

IF(  I. EG. MI    GO    TL    60 

AIMaAI [ ,M) 

DG     30    J=1,NN 
50    A(  I  ,J)=A (  I  ,J).-AIM*A{M,  J) 
60    CCN  riNUH 

GC     10    100 
7  0    DO    SO    I=M,N 

IF(A( I ,M ) .EG.  CO)    GC    TC    90 
C  DO    LOO?    80    EFFECTS    INTERCHANGE    OF    THf 

DO    SO    J  =  l,\i\ 

DUMY=A< I , J) 

A(  I  ,J)  =  M,\,J) 
CO    A(K»J)»DUMY 

GO  TO  30 
<-/u  CGNTINUE 

Wk  IT"  (3,  2) 

GO  TC  120 


ITH  AND  HTH  ROWS 
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ICO  CONTINUE 

TRANSFERING  INVERSE  OF  A  IN  THE  PLACE  bF  A 

DC  11C  1=1 iN 

DO  11C  J  =  l ,N 

JN-J+N 
110  A(  I  ,J)=A( I ,JN) 
120  RETLKN 

EN  13 


88 

SUBROUTINE  INTGKC 

c 

THIS  SUBROUTINE  NORMALIZES  THE  ITERATE, 

c 

FINOS  ITS  RALCIGHS  QUOTIENT  AND  FORMS 

c 

THE  INTIGRANO  T CR  THE  PROBLEM  OF 

c 

HEAT  TRANSFER  U.  AN  INCOMPRESSIBLE  FLUID 

c 

IN  LAM INAR  MOT  ION 

DIMENSION  F( IOC)  ,G( 100 >  ,H( IOC)  ,R{10CJ 

COMMON  M,N,IT,NC,NE,NP,MP1,NP1,EI 

CO'MKON  R  ( 3 ,8) ,  0(  8,16  J  t  C  ( 100)  ,E  { 10)  ,P(  100) 

1 

FORMAT (AS  16. 8) 

2 

FORMAT!/) 

c 

NORMALIZATION  OF  IHE  ITERATE 

DO  20  1=1, NC 

IF(C( I ) .ME. 0.0)  GC  TO  3C 

2  0 

CONTINUE 

30 

CU=C( I ) 

DC  40  J=l,NC 

40 

C(  J)=CU)/CIJ 

WRITE (3,1)  (C(  1  )  ,  1  =  1, NO 

WRITE( 3,2) 

:       c 

CALCULATION  OF  RALEIGHS  QUOTIENT 

CALL  S  EM UL  (  C  ,  N C  ,  E  ,  N L  ,  G  ,  i\G  ) 

WR11E(3,'2) 

CALL  SEMUL  (  C  ,  N  C  ,  G  ,  NG  ,H  ,  IVH ) 

WK  I  Tc(3,2) 

F(1)=0. 

CALL  tNTGRA(H,NH,F,NF, 1 ) 

W  R  I T  E I  3 ,  2 ) 

DN  =  C. 

DC  SO  1=1, NF 

9  0 

DN  =  CN*F(  I  ) 

CALL  DIFFCR(C,NC»F,iSF,  1) 

WRITE (3, 2) 

CALL  DIFF£R(F,NF,h,NK,  1  ) 

KR I  ft  I  3, 2) 

dc  lc  1*1, m 

• 
■ 

90 

c 
c 
c 

SUGHCUTINE  BOUNU 

THr  SUBROUTINE  CCMJlTLiS  THE  F  VECTOR 

FC«  ALL  EOUILITKIUM  PROBLEMS  AND  FO.'<  THE 

FIRST  i-iCDC  OF  ALL  EIGENVALUE  PROBLEMS 

DIMENSION  Y( LOl 

COMMON  M , M , I T , NC »  t\£ / NP , ^ PI ,NP1 , £1 

COMMON  B IS ,8 ) , 0 1 e ,16) ,C ( LOO) ,£ ( 10) ,P( 100) 

1  FGRKAT(4E16.8) 

2  FGi;!1vAT(//?OX,ciHL:CLMC/) 

DG  10  I=1,M 
10  C(  I  )»BU,MP1) 

CQN«*l. 

DG  JO  I*lfM 

Y(  I  )  =  0. 

A  I  ^  1 
'  DG  20  JM*MP1,NC 

AJMaJM 

Y(  I  )=Y{  I  )+CC'M*C  (  J.v) 
20  CCN=CON*  (  -i  J*/  (  AJM+1  .-A  I ) ) 

C0N=1. 

DC  30  J=l,  I 

AJM=M+1-J 
30  CON«C0««AJM 

DC  40  I*NP1»M 

DC  AC  J»1,M 
AG  CI  I  )  =  C( I Hft< If J)*Y< J) 

WRITE! 3,2) 

WRITE (3,1)  (C(I  )t  1  =  1  fM) 

RETURN 

• 

• 

■ 

• 

ENL, 

. 
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C 

c 
c 


SUBROUTINE  MATKUL  (A','N,P) 

THIS  SUftRCUTltfc  MULTIPLES  A  SQUARE  MATRIX  A 

OF  CRUER  N  BY  A  GOLIMN  VECTOR  B.  THE  ANSWER 

APPEARS  IN  THE  COLUMN  R. 

D  I  MENS  ION  A  (  8 1  16  J  •  B-l  8  )  i  C  (  8  ) 
1  FCRMATCE16.8) 
3  FORMAT 120X,6HMATMUL) 

HRITE(3f3J 

DC  IOC  1  =  1,  N 
D(  I  )  =  C. 
DC  ICO  J  =  1,N 
ICG  D(  I  )=C(  I )+A(I, J)»D( J) 

DC  20C  i  =  i,r-; 

2CO  B(  I  )  =  C(  I ) 

WRI  IE(3f  1)  (B(  I  )  i  1  =  1  iN) 

RETURN 

ENC 
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SUBROUTINE    DIFFuH  IC  f  NC  ,C  ,ND,N) 
NTH    DERIVATIVE    CF    SERIES   C    =    SERIES    0 
DIMENSION   C(  IOC  )  f  C  (  J.00  ) 
1    FCRfcAT{4E16.8) 

3    FORMAT    (/20X,     IC    HDEKIVMIVE       /) 
Wi<  IT£(3|31 
AF»1. 

OC    10    1=2, N 
A  I  --  I 
10    AF=AF#AI 

IFtN.EQ.ll      AF«l. 
DO    20    I=1,UC 
IN*  I+N 
AI«I 
A  I  N  =  I  .\ 

D(  I  )=C(  IN) »AF 
20   AF.*AF«AIN/AI 
ND=NC-N 

WRITEOil)     (Dm,!*    l.NC) 
RETURN 
END 
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SUBROUTINE    INTGRA    tC»NC-,0,NO,N  3 

NTH    INMGRM.    OF    SERIES    C    =    S  E  * ;  I  L:  S    0 

DIMENSION    C( 100)  ,C( ICO) 
1    FORMAT (AE16. 8] 
3   FORMAT    (/20X,     11    KIN  riGRATIO.N      /) 

WRITE ( it  J) 

Nl=.\+l 

AF»1. 

DC     10     I=2,N 

A  I  =  I 
10    AF=AF*AI 

IF(N.EO.l)    AF»li 

NCeNC+N 

DO    20    1=1,ND 

IN=i+N 

AI=I 

A  IN* IN 

Dl  IN)=CI  [)/AF 
20    AF»AF*AIN/AI 

Wr!IT£(3,l)     (D(  I  )  ,  l  =  .NL,NfJ) 

RETURN 

END 
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C 
C 


SUBROUTINE    RESULT 

THIS    SUBROUTINE    GIVL:5     fHE    RESULTS    OF    AN 

EIGENVALUE      PROBLEM 

Cdvi'  ON    M»N , IT,  NC  ,  .\E  ,  NP  ,  NP  I  .  NP1  ,  E  I 

CpKttQN    l>.  (  R  ,  8  )  ,  0  I  P  ,  16  )  ,  C  (  100  )  ,  F  (  10  )  ,  P  (  100  ) 

1  F  C  *■'  iv  A  T  (  A  E  1  6  .  B  ) 

2  FORNATi/l 
HRiTEO'tl)    EI 
w R I T £ ( 2  .  1 )    EI 
WRITE (3,2) 

00  20     1=1, NC 

J 

IF (CI  I  )  .NC.0.0)  GC  TO  AC 
20  CONTINUE 
40  EI*C( I ) 

DC  60  J=I,NC 
6  0  C( J  )='".( J) /EI 

WRITE! 3,1)  (C(  I  )  »  1-1  »NC) 

WRI1E(2,1)  ICU),I"ltNC) 

RETURN 

END 
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An  iteration  technique  for  solving  two-point  boundary  pro- 
blems was  studied.   The  technique  yields  an  analytical  solution 
in  the  form  of  a  power  series  in  the  independent  variable.   This 
type  of  technique  had  been  previously  applied  to  the  problem  of 
buckling  of  columns  by  L.  Vianello  and  to  the  problem  of  critical 
speeds  of  rotating  shafts  by  A.  Stodola.   Here  the  method  was 
applied  to  a  larger  class  of  problems. 

A  variety  of  problems  of  both  equilibrium  and  eigen-value 
type  were  selected  for  the  study.   The  problems  varied  in  the 
order  of  complexity  of  the  differential  equations  governing  them. 
Some  of  them  were  simple  enough  to  have  solutions  in  closed  form. 
These  solutions  were  expanded  in  power  series  so  as  to  make  the 
comparison  direct,  and  therefore,  easy.   Others  did  not  have 
closed  form  solutions  so  that  it  was  necessary  to  solve  them  by 
some  other  method.   The  solutions  thus  obtained  were  used  to  form 
a  basis  for  comparison. 

In  the  case  of  eigenvalue  problems  the  method  led  to  the 
mode  corresponding  to  the  lowest  eigenvalue.  The  process  was 
then  modified  to  extract  modes  corresponding  to  higher  eigen- 
values.  The  orthogonality  condition  was  used  for  this  purpose. 

Lastly,  a  problem  of  practical  importance  was  selected.   The 
problem  was  typical  of  its  class.   The  governing  equation  had  a 
lower  order  derivative  term  (independent  of  eigenvalue)  present 
with  the  highest  order  derivative  term.   The  terms  had  non- 
constant  coefficients.   To  the  author's  knowledge,  a  closed  form 


solution  had  not  been  obtained  for  this  problem  at  that  time. 
However,  a  numerical  solution  with  very  good  accuracy  was  avail- 
able.  This  solution  was  used  for  comparison  with  that  obtained 
by  the  present  method. 


